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Abstract. We study the bipartite entanglement per bond to determine characteristic 
features of the phase diagram of various quantum spin models in different spatial 
dimensions. The bipartite entanglement is obtained from a tensor network 
representation of the ground state wave-function. Three spin-1/2 models (Ising, XY, 
XXZ, all in a transverse field) are investigated. Imaginary-time evolution (TEBD in 
ID, ‘simple update’ in 2D and 3D) is used to determine the ground states of these 
models. The phase structure of the models is discussed for all three dimensions. 
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1. Introduction 

It is most intriguing to observe that matter comes in different phases. Usually, by 
a change of temperature one may switch between the different phases. Close to zero 
temperature quantum fluctuations dominate over thermal fluctuations, and quantum 
systems undergo quantum phase transitions (QPTs) as a result of a change of an external 
control parameter [1]. 

Discovery of new and exotic quantum phases and critical points, such as 
topologically ordered phases [2], quantum spin liquids [3], deconfined quantum critical 
points [4, 5], which are not encompassed by the Landau’s theory of phase transitions, 
recently stirred significant interest into the study of QPTs. Moreover, it was found 
that from quantum information theory perspective quantum phases and QPTs can be 
distinguished and characterized in terms of quantum entanglement [6, 7]. Entanglement 
measures are able to determine critical properties of the systems, in particular the 
positions of the critical points [8]. 

In the present paper we study how a specific entanglement quantifier, the bipartite 
entanglement per bond [9], can be used for a fast and computationally rather inexpensive 
determination of quantum critical points and other characteristic features of the phase 
structure for various one-, two- and three-dimensional quantum spin models. This 
quantifier does not require the calculation of expectation values from the states but 
is obtained directly from the representation of a state. This is the reason why this 
quantifier comes at low numerical cost. Calculation of expectation values for two 
and three-dimensional systems is numerically expensive, that is why this quantifier 
is particularly useful for the analysis in 2D and 3D. The spin-1/2 3D XY model in a 
transverse field is studied systematically here for the first time. 

The entanglement per bond is obtained directly from a tensor network 
representation of the ground state wave-function. Tensor networks (TN) [10, 11] provide 
a modern and promising tool for the numerical investigation of many-body systems. The 
basic idea of TN methods is to represent the wave function of a many body quantum 
system by a network of interconnected tensors. For a recent review see [12]. For one¬ 
dimensional (ID) systems the best known tensor network states are matrix product 
states (MPS) [11]. Their direct generalization to 2D and 3D are the projected entangled 
pair states (PEPS) [13]. Details about MPS and PEPS can be found in [14-17]. 

For our investigation we choose three spin-1/2 quantum models with nearest 
neighbor interactions and analyze them in the thermodynamic limit: the XY model 
(with quantum Ising model as a special case) and the XXZ model both in a transverse 
magnetic field. The models are investigated in one- (spin ring), two- (square lattice), 
and three (cubic lattice) dimensions. We use the imaginary-time evolution for the 
determination of the ground states in the TN representation. The entanglement per 
bond is calculated directly from the MPS and PEPS representation of the ground states. 

The paper is organized as follows. In Sec. 2 we review the definition of the 
entanglement per bond and comment on its physical meaning. In Sec. 3 we present 
numerical results along with a detailed discussion. Conclusions are made in Sec. 4. 
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2. Bipartite entanglement per bond 

In this section we define the bipartite entanglement per bond S'pb in terms of the bond 
vectors characterizing the TN. We comment on its physical interpretation and its ability 
to describe the entanglement properties of the state. Strictly speaking S'pb is not an 
experimentally measurable quantity, since it depends on a specific representation of a 
quantum state as a TN. 

Ground states in TN representation are obtained either variationally or via 
imaginary-time evolution [18] based on a Trotter expansion [19] of the evolution 
operator. In the present paper we use time-evolving block decimation (TEBD) [20] to 
determine ground states of one-dimensional systems and the ‘simple update’ scheme [16] 
in 2D and 3D. The algorithms we use are well-known, so we do not provide details of 
their implementation. 

The TEBD algorithm naturally leads to a ground state in the canonical form of 
MBS [14, 21, 22]. In this form, besides tensors A at each site one also has bond vectors 
A at each bond of the tensor network. Within the given canonical MPS the bond states 
(bond vectors) can be regarded as renormalized bases of the physical degrees of freedom 
of the many body system (e.g. ‘effective’ spins). Thus bond vectors can be treated as 
Schmidt coefficients in the wave function decomposition, i.e. they are objects which 
describe entanglement within the state. Moreover, the reduced density matrix of the 
state can be expressed approximately via virtual (bond) degrees of freedom, i.e. via 
the bond vectors. Unfortunately, no exact canonical form for PEPS exists, but the 
‘simple update’ [16, 23] algorithm for imaginary-time evolution leads to a ground state 
in approximate canonical form of PEPS. Thus, bond vectors can be used to obtain the 
bipartite entanglement per bond in this case as well. 

The numerical cost of the TEBD algorithm is 0{nA) with m the bond dimension. 
The numerical costs for ‘simple update’ procedure in 2D and 3D are 0{m^) and 0{m^^) 
respectively. Notice, that the calculation of expectation values from a state given in MPS 
or PEPS representation requires numerically expensive (especially in 2D and 3D) tensor 
contraction procedure, but the bipartite entanglement per bond is obtained directly 
from the wave function tensor network representation, which makes this quantifier so 
attractive. 

The bipartite entanglement per bond S'pb is directly defined in terms of the 
components of the bond vectors A, which are normalized to satisfy A^ = 1. Then, 
S'pb is defined as the entanglement entropy [24] or von Neumann entropy 

•S'pb = (1) 

i 

In practice, we find that the bond vectors connecting to a given site i are approximately 
equal, which is a consequence of the translational symmetry of the tensor network we 
implement. As a consequence, we average the bond vectors connecting to a site i for 
the calculation of S'pb. 

Note that the maximally possible value for S'pb (measuring the maximally possible 
entanglement in the state) is dependent on the chosen bond dimension m of the MPS 
or PEPS representation. This means that such a representation may not be able to 
accurately describe states close to critical points where entanglement may be very 
large. Still, the calculated entanglement per bond provides at least semi-quantitative 
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information on the phase strncture of the many-body system. 

3. Numerical results and interpretation 

The present section provides nnmerical results for the bipartite entanglement per bond 
for the Ising, XY and XXZ models in a transverse held in one, two, and three dimensions. 
Ground states are obtained via imaginary-time evolution algorithms: time-evolving 
block decimation (TEBD) for modeling ground states of one-dimensional systems and 
the ‘simple update’ scheme for two- and three-dimensional systems. We choose the spin 
ring for one-dimensional geometry, the square lattice in 2D and the cubic lattice in 3D 
always with periodic boundary conditions and translationally invariant. Critical points 
obtained from this measure dependence are compared to results from previous studies, 
based on other algorithms. 

Imaginary-time steps are proceeded until the convergence of bond vectors is reached. 
As a check for a good ground state we use the condition of approximate equality of all 
bond vectors, that means presence of rotational symmetry in the ground state. Time 
step size in the imaginary-time evolution is reduced repeatedly during convergence. For 
TEBD we use MBS bond size up to m = 40. For the ‘simple update’ algorithm we take 
bond sizes m = 4 and m = 2 in 2D and 3D, correspondingly. 

3.1. Anisotropic XY model 

We hrst study the anisotropic XY model in a transverse magnetic held, 

- E {(1 + 7)Sf» sj + (1 - US’ 0 + hY.s‘, (2) 

hd) * 

where {i,j) indicates a summation over nearest neighbors. It includes the Ising model 
as a special case (7 = 1). We use periodic boundary conditions, and the spin operators 
S'" = |cr" are related to the Pauli matrices a". The model has two parameters, the 
anisotropy 7 and the magnetic held h. 

The one-dimensional XY model can be solved analytically [1, 25, 26]), and as guide 
for the reader its well-known phase diagram is shown schematically in Figure 1. There 
are ferromagnetic and paramagnetic phases which are separated by the black dashed 
critical lines. On the red circle -|- 7 ^ = 1 the model is classical, where we expect 
entanglement to vanish. It is known that correlation functions have an oscillatory 
tail [27] inside this circle. Consequently, this region often is called ‘oscillatory’. The line 
7 = 0 separates two diherent ferromagnetic phases: the x-phase (upper half plain) and 
the y-phase (lower half plain). 

The phase structure of the XY model in 2D is similar to the one shown for ID [9, 28]. 
It is classical for (h/2)^ -|- 7 ^ = 1. For more details on the phase diagram of ID and 2D 
anisotropic XY models see [28, 29]. 

The 3D spin-1/2 anisotropic XY model in a transverse held has not yet been studied 
in detail according to our knowledge. Our results suggest that the 3D XY model has 
a similar phase structure as its ID and 2D counterparts, and that it is classical on the 
circle (/i/3)^ -I- 7 ^ = 1. However, the bipartite entanglement per bond does not provide 
enough information to determine the detailed characteristics of the phases. In order to 
do so, one has to calculate expectation values of various observables. The entanglement 
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Figure 1. Phase diagram of the ID XY model. Black dashed lines indicate phase 
separation lines. Green dashed lines 7 = 1, 7 = 0.7, 7 = 0.5 and h = 0.5 indicate the 
cases presented in Figs. 2, 3, 4, and 5. The dots mark points where S'pb shows special 
characteristics. (PM: Paramagnetic phase, FM: Ferromagnetic phase). 


per bond S'pb just provides information about the position of critical points as well as 
lines of vanishing entanglement (‘classical’ lines). 

We now discuss results for S'pb along various lines in the phase diagram (indicated 
as green dashed in Fig. 1): 7 = 1 (Ising model), 7 = 0.7 and 7 = 0.5 as well as h = 0.5. 

The results for 7 = 1 are shown in Figure 2. In ID our TEBD calculations provide 
a critical point in very good agreement with the analytical result: hl^ ^ 1.00. The 2D 
quantum Ising model cannot be solved analytically, and various methods are applied 
to solve it numerically, e.g. quantum Monte-Carlo (QMC) methods. Such calculations 
hnd a transition between a ferromagnetic and a paramagnetic phase at a critical point 
hf^ = 3.044 [30]. The tensor network implementation applied here produces numerical 
results significantly faster than QMC calculations, however, with less precision; it 
determines a critical point at hf^ ~ 3.26. More precise results can be obtained with 
more elaborate tensor network implementations and larger bond sizes [12]. However, 
it is our goal to investigate S'pb properties using small numerical cost, i.e. with small 
virtual bond dimension. 

The 3D quantum Ising model cannot be solved analytically as well. The series 
expansion study in the T —)■ 0 limit [31] predicted the critical point of hf^ ^ 5.14. 
Another study based on ‘simple update’ scheme was performed in [17] and the critical 
point was determined at h^^ ^ 5.29 from the ground state magnetization. The 
magnetization calculations presented in [17] require calculation of expectation values, 
that is the tensor network contraction procedure. We obtain a result {hf^ ^ 5.30) in 
good agreement with latter one, but without the need to calculate a matrix element 
explicitly. 

Numerical results for 7 = 0.7 and 7 = 0.5 are shown in the Figures 3 and 4, 
respectively. S'pb has two characteristic features: at the ferromagnetic-to-paramagnetic 
phase transition it shows a maximum, and it is zero at the boundary of the oscillatory 
region {h/d^ -|- 7 ^ = 1 (d = 1, 2, 3). Both points are easily identified in the hgures. 
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Figure 2. Bipartite entanglement per bond for 7 = 1 (Ising model) in ID (chain), 
2D (square lattice), 3D (cubic lattice) as a function of the magnetic field h. 2D and 
3D results are multiplied by a factor of 10 and 15, respectively. Bond dimensions: 
= 20, = 4, = 2. 


In the case 7 = 0.7 (Figure 3) characteristic features are seen at hl^ ~ 0.71 and 
at hl^ ~ 1.00 for ID. These results perfectly agree with theoretical predictions [9, 29]. 
The 2D result 1.42 perfectly agrees with the prediction from (h/d)^ + 7 ^ = 1, 

while h® ~ 2.89 differs slightly from the reference value h® = 2.72 obtained from finite 
size scaling [28]. The 3D boundary point hf^ ~ 2.14 satisfies (h/3)^ + 7 ^ = 1, and we 
obtain the 3D ferromagnetic-to-paramagnetic critical point at hf^ ~ 4.63. 



0 1. 2 2.89 4 4.63 5 6 


h 

Figure 3. Bipartite entanglement per bond for 7 = 0.7 in ID (chain), 2D (square 
lattice), 3D (cubic lattice) as a function of magnetic field h. 2D and 3D results 
are multiplied by a factor of 12 and 24, respectively. Bond dimensions: = 20, 

= 4, = 2. 

For 7 = 0.5 we obtain in ID h^^ ~ 0.87 and hl^ ^ 1.00 again in perfect agreement 
with theoretical predictions. Our 2D result h^^ ^ 1.73 agrees with the theory while 
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hfP K, 2.64 differs slightly (the reference value is he = 2.5 [28]). We observe zero 
entanglement at hfP ~ 2.60, which satishes the expression (/i/3)^ + 7 ^ = 1. The 
3D ferromagnetic-to-paramagnetic critical point for 7 = 0.5 is determined to be at 
hfP ^ 4.17. 



0 1 2 2.64 4.17 5 6 


h 

Figure 4. Bipartite entanglement per bond for 7 = 0.5 in ID (chain), 2D (square 
lattice), 3D (cubic lattice) as a function of magnetic field h. 2D and 3D results 
are multiplied by a factor of 12 and 24, respectively. Bond dimensions: = 20, 

to^d _ _ 2 

Figure 5 represents bipartite entanglement per bond for the XY model in ah three 
dimensions as a function of anisotropy parameter 7 at hxed magnetic held h = 0.5. 
The entanglement measure peaks at 7 = 0 indicating the phase transition between two 
different ferromagnetic phases and vanishes at points which correspond the equation 
(/i/d)^ + 7 q = 1. Figure 6 displays an enlarged view of the region where entanglement 
vanishes for positive 7 . The obtained values for 7 ^° ~ ±0.866, 7 ^° ~ ±0.968, 
780 ^ ±0.986 agree to a high precision with the expected theoretical values. 

The boundary points of the oscillatory region are obtained numerically up to a 
very high precision as predicted by the relation (/i/d)^ ± 7 ^ = 1. This observation nicely 
underlines that MPS and PEPS tensor network states are particularly able to model 
ground states with a low amount of entanglement, that is to model states for which a 
low bond dimension is enough. 

On the other hand, our results for the phase boundaries which are characterized 
by a maximum value of Spb differ from the theoretical predictions in 2D presented in 
Ref. [28]. This is due to the low virtual bond dimension of the PEPS used here, i.e., 
its impossibility to capture the large amount of entanglement needed to determine the 
critical point more precisely. Consequently, we would expect that the results for the 3D 
phase boundaries are even less precise than our 2D results. 

From the Figures shown above in this section we conclude that S'pb is quite capable 
of determining important features of the phase diagram for a wide range of model 
parameters at relatively small numerical cost. The results for the ferromagnetic-to- 
paramagnetic phase transition in the 3D XY model, namely hf^ ~ 4.17 for 7 = 0.5 and 
^3D ^ ^ ^ _ g Were not obtained before to the best of our knowledge. 












Phase diagram of one-, two-, and three-dimensional quantum spin systems 



-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 

r 

Figure 5. Bipartite entanglement per bond for h = 0.5 in ID (chain), 2D (square 
lattice), 3D (cubic lattice) as a function of anisotropy parameter 7 . 2D and 3D results 
are multiplied by a factors 18 and 30, respectively. Bond dimensions: m}^ = 40, 
= 4, 771^° = 2. 
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Figure 6. Bipartite entanglement per bond for h=0.5 in ID (chain), 2D (square 
lattice), 3D (cubic lattice) as a function of anisotropy parameter 7 . Displayed region 
shows those positive values of 7 , where entanglement vanishes. ID, 2D and 3D results 
are multiplied by a factors 15, 200, and 400, respectively. Bond dimensions: = 40, 

= 4, = 2. 


3.2. XXZ model 

Next we study the spin-i XXZ (anisotropic Heisenberg) model in a transverse magnetic 
field, 

hXxz = ^ 1^- (g, + Sl®S] + XSt® S';} -Y.hSt (3) 

(hi) * 

as a function of the anisotropy parameter A and the magnetic field h. 

The phase structure of the XXZ model was studied extensively in ID [32-37], 
2D [35, 38, 39], and 3D [38]. In all three dimensions the model shows three phases [38]: 
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an antiferromagnetic (Neel) phase, a XY (spin-flopping) phase and a ferromagnetic 
phase. These three phases are separated by two critical lines, and we denote them as he 
and hg. For ID the well-known phase diagram can be obtained analytically using the 
Bethe Ansatz [32-34], and it is shown in Figure 7. The ferromagnetic phase is separated 
from spin-flopping phase by the line hg = d {1 -\- A) (d = 1, 2, 3 is the dimension of the 
model) [38, 40, 41], In ID the Neel phase is separated from the spin-flopping phase by 
the curve he [34, 42] 


he 


TT sinh A 
A 


E sech|^(l + 2n) 


(4) 


with A = arccoshA. In particular, at zero magnetic field the ferromagnetic-to-XY 
critical point is at A = —1 and the XY-to-antiferromagnetic critical point at A = 1. 

Similar phase diagrams for the two- and three-dimensional XXZ models, based on 
quantum Monte Carlo studies, can be found in [38]. The most notable difference to the 
ID phase diagram shown here is the fact that the critical lines he have a finite derivative 
with respect to h at A = 1. 



A 

Figure 7. Phase diagram of the ID XXZ model. The line hg separates the XY phase 
from the ferromagnetic phase (FM). The curve he separates the anti-ferromagnetic 
(AFM) from the XY phase. The green dashed lines h = 1 and A = 1.5 correspond 
to the cases discussed in detail in the text and in Figures 9 and 10. The black dots 
indicate the critical points obtained from the entanglement per bond. 


For our numerical studies we choose three cases: h = 0 (Figure 8) and h = 1 
(Figure 9) as well as A = 1.5 (Figure 10). These lines are denoted as green dashed lines 
in the phase diagram 7. 

Figure 8 presents our results for the bipartite entanglement per bond at zero 
magnetic field as a function of A. From the figure we see that critical points A = ±1 are 
obtained correctly. Notice that the critical point A = 1 is characterized by a maximum 
of the entanglement, but that at the critical point A = — 1 the entanglement vanishes. 
Moreover we observe that at A = —1 App shows a discontinuity in ID, while it increases 
slowly in 2D and 3D. 

Figure 9 shows our results for the bipartite entanglement per bond as a function of A 
at h = 1. The ID critical points are in a very good agreement with exact values [33, 34]: 
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A 


Figure 8 . Bipartite entanglement per bond for the XXZ model with zero magnetic 
field in ID (chain), 2D (square lattice), 3D (cubic lattice) as a function of anisotropy 
parameter A. 2D and 3D results are multiplied a factor of 5 and 6, respectively. Bond 
dimensions: = 20, = 4, = 2. 


the ferromagnetic to spin-flopping phase transition occurs at = 0 and the spin- 
flopping to antiferromagnetic transition at A^° = 2.75. The 2D critical points are 
found at A^^ = —0.5 and A^° = 1.18, and the 3D critical points at A^^ = —0.67 
and A^° = 1.07. The values obtained for A^ in all three dimensions correspond to the 
relation hs = d (1 -|- A) with dimensions d = 1, 2, 3. Our values obtained for Ac for the 
two- and three-dimensional XXZ model are in a good agreement with quantum Monte 
Carlo results given in [38]. 



A 


Figure 9. Bipartite entanglement per bond for XXZ model with h = 1 in ID 
(chain), 2D (square lattice), 3D (cubic lattice) geometries as a function of anisotropy 
parameter A. 2D and 3D results are multiplied a factor of 8 and 12, respectively. Bond 
dimensions: = 20, = 4, = 2. 
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Figure 10 shows results for Spb as a function of the magnetic held h at A = 1.5. 
The XXZ model at A = 1.5 was studied analytically in ID [33, 34] and was investigated 
numerically in 2D and 3D [9, 38]. These studies predict two critical points for each 
dimensionality. From equation hg = d (1 + A) one obtains = 2.5, hf^ = 5.0, 
= 7.5, and our results agree with these values. For the one-dimensional model the 
exact solution gives hfP ~ 0.0866 ([33]), and our result hfP ~ 0.09 is very close. Our 
2D and 3D results are hfP k. 1.8 and hfP ra 3.2 are in a good agreement with previous 
quantum Monte Carlo studies [38]. 

From Figures 8, 9, and 10 we conclude that Spb is capable of providing quite some 
insight into the phase structure of a quantum spin model. One not only obtains critical 
points, but also points or lines where the system behaves classically. The variation of 
SpB in between such lines provides insight how the entanglement of the various ground 
states changes as a function of the control parameters. Of course, for a more detailed 
investigation one then would need to calculate relevant expectation values in order to 
obtain physical properties in interesting regions. The identihcation of such regions is 
easily guided by the Spb, which is obtained at a rather low numerical cost. 



Figure 10. Bipartite entanglement per bond for the XXZ model at A = 1.5 in ID 
(chain), 2D (square lattice), 3D (cubic lattice) as a function of external magnetic field 
h. 2D and 3D results are multiplied by factors 7 and 12, respectively. Bond dimensions: 

= 20, = 4, rrP^ = 2. 

From Figures 8, 9, and 10 we observe interesting differences for S'pb at the critical 
points in different spatial dimensions. E.g., at the transition from the ferromagnetic 
to the spin-flopping phase S'pb shows a jump in ID, for 2D and 3D it shows a cusp. 
The numerical results indicate that the entanglement structure of the ground state of 
the XXZ model in different spatial dimensions differs not only quantitatively bnt also 
qnalitatively. Understanding the reasons for these differences requires further analysis 
which is beyond the bipartite entanglement per bond quantiher. 
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4. Conclusions 

In the present paper the bipartite entanglement per bond S'pb was calculated for the 
Ising, anisotropic XY and XXZ models in a transverse magnetic held and in one, two, 
and three dimensions for various sets of control parameters. All obtained results are in 
very good agreement with previous studies. They show that bipartite entanglement per 
bond is an efficient, easy-to-calculate and relatively precise tool for the determination 
of the some characteristics of the phase diagram of quantum models. It conveniently 
uses the tensor network representation of a state. 

The three-dimensional anisotropic XY model with transverse held on a cubic lattice 
was studied numerically for the hrst time. The obtained boundary points of the 
oscillatory region satisfy the relation (ho/3)^ -|- 7 ^ = 1 . 

We observe that the bipartite entanglement per bond determines characteristic 
features of the phase diagram better if states with a low amount of entanglement are 
involved. For our 2D and 3D calculations we took rather low virtual bond sizes and 
therefore our PEPS can only model states with rather low entanglement. Moreover, 
if more detailed questions about the phase diagram need to be answered for specihc 
regions (e.g. the magnitude of the magnetization), then calculations of matrix elements 
cannot be avoided. 

The present paper was geared towards a semi-quantitative analysis of the bipartite 
entanglement per bond. An interesting and promising task is the quantitative 
investigation which requires implementation of the symmetries of the ground states 
into the imaginary-time evolution algorithms. That would enable higher virtual bond 
dimensions for MPS and PEPS and the calculation of the entanglement spectrum of the 
state. 

Besides the geometries studied in the present paper (spin ring, square lattice, 
cubic lattice), the bipartite entanglement can also be used to characterize properties 
of quantum models on hyperbolic lattices [43], after modihcation of the ‘simple update’ 
algorithm for the ground-state calculation. Moreover, the present study should be 
extended to spin -1 systems. 
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